The Extended Bose-Hubbard Model on a Honeycomb Lattice 
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We study the extended Bose-Hubbard model on a two-dimensional honeycomb lattice by using 
large scale quantum Monte Carlo simulations. We present the ground state phase diagrams for both 
the hard-core case and the soft-core case. For the hard-core case, the transition between p = 1/2 
solid and the superfluid is first order and the supersolid state is unstable towards phase separation. 
For the soft-core case, due to the presence of the multiple occupation, a stable particle induced 
supersolid ( SS-p ) phase emerges when 1/2 < p < 1. The transition from the solid at p = 1/2 
to the SS-p is second order with the superfluid density scaling as p s ~ p — 1/2. The SS-p has the 
same diagonal order as the solid at p — 1/2. As the chemical potential increasing further, the SS-p 
will turn into a solid where two bosons occupying each site of a sublattice through a first order 
transition. We also calculate the critical exponents of the transition between p — 1/2 solid and 
superfluid at the Heisenberg point for the hard core case. We find the dynamical critical exponent 
z = 0.15, which is smaller than results obtained on smaller lattices. This indicates that z approaches 
zero in the thermodynamic limit, so the transition is also first order even at the Heisenberg point. 



Bose-Hubbard model and its derivatives have been ex- 
tensively studied as the low energy effective models for 
ultracold atoms in optical lattices where the superfluid 
to the Mott insulator transition exists in the zero tem- 
perature 0, 0, Q , as well as models for the possible su- 
persolid phases in both bipartite and frustrated lattices 
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the Fermi Hubbard model, the Bose-Hubbard model can 
be used to explore the effects of the strong correlations. 
Unlike the fermi system, there is no "sign problem" in 
bosonic case, so the Quantum Monte Carlo (QMC) sim- 
ulations can be successfully performed. Many interesting 
results have been found in these models, such as super- 
solid phase, valence-bond-solid phase and striped phase. 

In the past several years, although the extended Bose- 
Hubbard model on two-dimensional square lattice [H, 0] 
and triangular lattice [7|, |8|, |lfj, LUj has been studied 
extensively using QMC methods, another common lat- 
tice, honeycomb lattice has rarely been studied by QMC. 
In [l3|, by using a dual vortex method (DVM), one of 
the authors studied phases such as superfluid, solid, su- 
persolid and quantum phase transitions in an extended 
boson Hubbard model such as Eqn[T] slightly away from 
half filling on bipartite lattices such as honeycomb and 
square lattice. He found that in the insulating side, differ- 
ent kinds of supersolids are generic possible states slightly 
away from half filling and also proposed a new kind of su- 
persolid: valence bond supersolid which maybe stabilized 
by possible ring exchange interactions. He showed that 
the quantum phase transitions from solids to supersolids 
driven by a chemical potential are in the same universal- 
ity class as that from a Mott insulator to a superfluid, 
therefore have exact exponents z = 2,u = 1/2, -q = 0. 
For example, near the solid to the supersolid transition, 
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FIG. 1: The ground-state phase diagram for 2D hard-core 
Bose-Hubbard model on the honeycomb lattice in the grand 
canonical ensemble obtained from quantum Monte carlo sim- 
ulation. SF: superfluid phase. The solid line is first order 
transition, while the dashed line is second order. 



the superfluid density inside the supersolid phase should 
scale as p s ~ \p-l/2\( d+z - 2 > = \p-l/2\ = 5f with possi- 
ble logarithmic corrections. Comparisons with previous 
quantum Monte-Carlo (QMC) simulations on a square 
lattice are made. 

The DVM is a magnetic space group ( MSG ) 
symmetry-based approach which can be used to classify 
some important phases and phase transitions. But the 
question if a particular phase will appear or not as a 
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FIG. 2: (color online), (a) The boson density p as a function 
of chemical potential p for /? = 100V and L = 24. (b) and 
(c) are superfiuid density p s and pa-b as functions of p for 
t/V = 0.4, = 100V, and L = 12, 24, 30 respectively. 



ground state can not be addressed in this approach, be- 
cause it depends on the specific values of all the possible 
parameters in the EBHM in EqnfTJ So a microscopic ap- 
proach such as Quantum Monte-Carlo (QMC) is needed 
to compare with the DVM. The DVM can guide the QMC 
to search for particular phases and phase transitions in a 
specific model. Finite size scalings in QMC can be used 
to confirm the universality class discovered by the dual 
approach. So far, there is no QMC simulations on a hon- 
eycomb lattice. In this paper, we study the extended 
Bose-Hubbard model on a honeycomb lattice by QMC 
simulations and also compare with the results achieved 
in [l3| by the DVM. We find that the ground-state phase 
diagram is qualitatively similar to the one obtained on 
square lattice. In the hard-core limit, the supersolid state 
is unstable towards the phase separation. The transition 
between p = 1/2 solid and the superfiuid is the first order 
one. In the soft-core case, due to the presence of the mul- 
tiple occupation, a stable particle doped supersolid ( SS-p 
) phase emerges when bosons are doped into p = 1/2 solid 
(i.e., fillings p > 1/2). We find the solid at p = 1/2 to the 
SS-p at p > 1/2 is a second order transition, the super- 
fluid density inside the SS-p scales as p s ~ p— 1/2 which 
is consistent with the result achieved in [13fl by the DVM. 
Very precise finite size scalings by QMC [T7J are under- 



way to test the exact exponents z = 2, v = 1/2, 77 = 0. 
However, a hole doped supersolid phase remains unsta- 
ble to phase separation when vacancies are doped into 
p = 1/2 solid (i.e., fillings p < 1/2). We also calculate 
the critical exponents of the transition between p = 1/2 
solid and superfiuid at the Heisenberg point. We find the 
dynamical critical exponent z and the correlation length 
exponent v are z — 0.15 and v = 0.38. The dynami- 
cal critical exponent we have obtained is smaller than it 
was previously obtained [IS]. We speculate that z will 



approach zero as the size of the simulated system goes 
to infinity, so the the SF to the solid transition at the 
Heisenberg point remains first order. 

The extended Bose-Hubbard model we study can be 
written as: 

h = -t ^2 ( a l a 3 + a ] a i) + 1* n ^ nt ~ 1 ) 



<»j> 

<ij> i 



(1) 



where aj(aj) is the creation (annihilation) operator of 
the bosonic atom at site i; m = a\a,i is the occupation 
number and p is the chemical potential. < ij > runs 
over nearest neighbors. U and V represent the on-site 
and nearest-neighbor repulsive interaction, respectively. 
When U/t — > 00, the Hamiltonian reduces to the hard- 
core one with no multiple occupation. The method we 
employ is the stochastic series expansion Quantum Monte 
Carlo (QMC) method with directed loop (SSE) 0. 

In order to characterize different phases, usually one 
studies the static structure factor S(Q) and superfiuid 
density p s 20], 
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where W is the winding number fluctuation of the 
bosonic lines, (3 is the inverse temperature, and N = 
2 x L x L is the lattice size. In a honeycomb lattice, 
S(Q = (47r/3, 0)) is always nonzero when p > 0, because 
Q = (47r/3,0) is not only the wave vector for the solid 
phase but also the reciprocal vector on a honeycomb lat- 
tice. So it can not be used to characterize the solid order 
on a honeycomb lattice. Here we use the density differ- 
ence between the two sublattice, pa-b, to characterize 
the solid order for a honeycomb lattice, 



PA-B 



\PA ~ PB\ 



(3) 



where pa and pb are the boson density for A and B 
sublattice respectively. The two parameters p s and pa-b 
can characterize different phases in the system: a solid 
phase is characterized by pa-b ^ and p s — 0; a SF 
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phase by pa-b = and p s ^ 0; a SS phase by both 
Pa-b 7^ and p s ^ 0. Different solid phases can be 
categorized by different values of pa-b- 

We first consider the hard-core limit, which is rela- 
tively simple. The phase diagram is presented in Fig. 
1. This phase diagram can be qualitatively explained by 
strong-coupling arguments. For p < —Zt, the system is 
empty while for p, > 3(i + V), it is a Mott insulator with 
one boson per site. At large values of t/V, the system 
exhibits a supcrfluid phase. A solid phase with p = 1/2 
emerges in small values of t/V. At p/V — 1.5, the crit- 
ical point between superfluid phase and p — 1/2 solid is 
at {t/V) c = 0.5. This point (t/V, p/V) = (0.5,1.5) is 
the Heisenberg point. 

When bosons (holes) are doped into the p = 1/2 solid, 
the domain wall proliferation mechanism, i.e., the addi- 
tional bosons (holes) can hop freely across the domain 
wall, gain the kinetic energy linearly in t, hence excludes 
the stable supersolid state [a, @] ■ The transition from the 
p = 1/2 solid to the superfluid is the first order one (see 
Fig. 2(a)). The phase diagram is qualitatively similar 
to the one obtained for a square lattice. This is not sur- 
prising since both the square lattice and the honeycomb 
lattice are bipartite lattices. However, the phase bound- 
ary is quantitatively different. This comes from the dif- 
ferent coordinate number of two lattices, i.e., q — 3 for 
the honeycomb lattice and q = 4 for the square lattice. 
Fig. 2 (a) shows the density p as a function of chemical 
potential p. For small values of t/V (say t/V = 0.2, 0.4), 
there is a jump from p = 1/2 to p > 1/2, which indicates 
a first-order transition. In the grand canonical ensemble, 
a jump in p is the token of a PS region. Fig. 2 (b) and 
(c) are the superfluid p s and pa-b as functions of p for 
t/V = 0.4. For p/V < 2, p s = while p A -B / 0, the 
ground-state is the p = 1/2 solid; for 2 < p/V < 4.2, 
p s ^ while pa-b — 0, the ground-state is a super- 
fluid phase; at the critical point pc/V = 2, both p s and 
PA-B have a jump, indicating the first order transition; 
for p/V > 4.2, both p s = and pa-b = with p = 1, 
the ground-state is a Mott insulator. There is no region 
for both p s and Pa-b nonzero, i.e., no stable supersolid 
state. 

In the above we have found that there is no stable 
supersolid phase in the hard-core Bose-Hubbard model 
in a honeycomb lattice, similar to a square lattice. It was 
found that by softening the onsite interaction U, i.e., in 
the soft-core case, a stable supersolid state emerges on 
the 2D square lattice [f|. Motivated by this, we next 
consider the soft-core case on the 2D honeycomb lattice. 
Hereafter we fix U = 15i and (3 = 2L in our calculations. 

Fig. 3 shows the phase diagram for the soft-core honey- 
comb lattice at U = 15i. The phase diagram is not sym- 
metric respect to p/V = 1.5 because finite U breaks the 
particle-hole symmetry. Besides the p = solid, p = 1/2 
solid (we named it CDW I here) and the Mott insulator 
(U > 3V) which already emerge in the hard-core case, an- 
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FIG. 3: The ground-state phase diagram for 2D soft-core 
Bose-Hubbard model on the honeycomb lattice obtained from 
quantum Monte carlo simulation. SF: superfluid phase; SS: 
supersolid phase; PS: phase separation; CDW I: p = 1/2 solid 
with one sublattice occupied and one boson per site; CDW II: 
p — 1 solid with one sublattice occupied and two bosons per 
site; MI: Mott insulator phase with one boson per site. The 
solid line is first order. The dashed line is second order. 



other p = 1 solid (U < 3V) with one sublattice occupied 
by two bosons per site emerges. We call this phase CDW 
II. A supersolid state emerges when atoms are doped into 
the p — 1/2 solid, while the supersolid phase is unstable 
towards the phase separation when holes are doped into 
the p = 1/2 solid due to the domain wall mechanism 
just like what happens in the hard-core case. The main 
reason for the emergence of the CDW II and supersolid 
state, is the possible multiple occupation in the soft-core 
case. When bosons are doped into the p = 1/2 solid, 
the additional bosons can occupy either occupied or un- 
occupied sites. When \3V — U\ ~ t, the doped bosons 
can form a superfluid on top of the p = 1/2 solid back- 
ground, this phase is a supersolid phase. Similar results 
have been found in the 2D soft-core square lattice , 2D 
triangular lattice ll| and ID chain 21 1 . 

Fig. 4 (a) is the plot of p v.s. p for U — 15i.For a 
density p just less than 1/2, there are jumps in p — p plot, 
which indicate the transition is a first order transition 
and the ground-state is a phase separation. This also 
indicates that the possible hole-doped supersolid ( SS-h 
) phase is unstable against the phase separation. The 
jump also emerges in the fillings just less than 1 for large 

V (say V = 6i). This indicates that the CDW II to the 
supersolid transition is 1st order. Fig. 4 (b) and (c) 
shows the finite size scaling of superfluid density p s and 
the diagonal long range order pa-b as functions of p for 

V = At and 6t, respectively. From these figures, we can 
identify the regions of superfluid with p s ^ 0; p = 1/2 
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FIG. 4: (color online), (a): The boson density p as a function 
of p for V = it, 5t, 6t at U = 15t and (3 — 2L for soft-core 
case, (b) (c): The finite size scaling of superfluid density p s 
and diagonal long range order pa-b as a function of p for 
V ~ it, 6i respectively for soft-core case. 



solid with pa-b ^ 0; supersolid with both p s ^ and 
Pa-b 7^ 0; phase separation with p s and pa-b showing 
discontinuity. In the inset of FigJUc) we also find that 
near p = 1/2, superfluid density scales as p s ~ (p — 1/2) 1 . 
This is consistent with the superfluid scaling derived by 
the DVM in [3]. pa-b also changes smoothly acrose the 
the CDW-I to the SS-p transition. So the SS-p has the 
same diagonal order as the CDW-I. These facts conclude 
that the transition from CDW I to SS-p is of second order. 
In fact, as shown by the DVM, the transition is in the 
same universality class as that from a Mott insulator to 
a superfluid, therefore has exact exponents z = 2, v = 
1/2,7] = 0. Very precise finite size scalings by QMC [ItJ 
are underway to test the universality class 17]. Possible 
ring exchange terms will also be added in 17f to test the 
possible valence bond supersolid proposed by the DVM 
in Q. 

As shown by the DVM in [13| , due to the change of sad- 
dle point structure of the dual gauge field on the CDW-I 
and SF, the SF to the CDW-I is a first order transition. 
The SF to the SS-p is also first order. Fig.l shows that 
the SF to the CDW-I is indeed first order. However, in 
Fig. 3, our QMC data in Fig. 3 can not determine if the 
SF to the SS-p is first order or second order. More pre- 
cise QMC such as double-peak histogram maybe needed 



to determine the nature of the SF to the SS-p transition. 

For the hard core case, at t/V = 0.5, corresponding to 
the isotropic Heisenberg point of the spin-1/2 quantum 
antifcrromagnetic model, the half filled solid melts to su- 
perfluid. The critical behavior of this Heisenberg point 
is controversial. By assuming a second order transition, 
Hebert et.al. [l8| found rather exotic critical exponents, 
z = 0.25 and v — 0.36. However, by studying the XXZ 
model crossing the isotropic point, Sandvik and Melko 
point out that the critical exponent z will become zero 
in the thermodynamic limit [22j . 

In the below, by taking the same strategy as Hebert 
et.al. [HI, e.g. assuming it is a second order quantum 
phase transition crossing the Heisenberg point, we also 
calculate the critical exponents for the 2D hard-core bo- 
son model on the honeycomb lattice at the Heisenberg 
point. The finite size scaling function of superfluid den- 
sity is assumed to have the following form 



Ps 



L 



l-F((t/V-(t/V)c)L^,0/L% 



(4) 



where z is the dynamical critical exponent, v is the corre- 
lation length exponent, and F is the corresponding scal- 
ing function. Note that the finite size scaling function F 
does not only depend on the scaled distance to the criti- 
cal point ((t/V — (t/V)c)L 1 / v ), but also on the ratio of 
the inverse temperature and the lattice size (/3/L z ). To 
obtain the critical exponents, one should simulate several 
sets of lattice sizes (L), each set with different values of z. 
Here we follow the procedure of Ref. [H[ for the square 
lattice, e.g. we fix a large enough (3 so that we can assume 
the second term of F is a constant as L is varied. With 
this choice, p s L z must be independent of L at the critical 
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FIG. 6: (color online). Data collapse of the superfluid den- 
sity p a for P = 100V at fi/V = 1.5. This yields the critical 
exponents z = 0.15 and v — 0.38. 

point (t/V = 0.5). We calculate (3 = 20,40,60,80, and 
100 and find the results depend little on temperature. In 
the below, we present results only of j3 = 100. 

First we simulate on smaller lattices L = 6,8, 10, 12 
and we obtain the same exponents as Ref. [3] (not show 
here). Then we simulate on larger lattices, e.g. L = 
16,18,24,30,36. The best intersection we found is at 
z = 0.15 (Fig. [5]). Replotting this figure by the scaled 
variable (t/V — 0.5)L 1 / l/ , the best data collapse we found 
is at v = 0.38 (see fig. [6]). 

The dynamical exponent we have obtained on larger 
lattices is smaller than that on smaller lattices. We spec- 
ulate that it will approach zero in the thermodynamic 
limit, so the transition is first order, though it maybe 
very weak. This scenario is consistent with Sandvik's 
conclusion on XXZ model on square lattice [22I ]. 

In conclusion, we have studied the extended Bose- 
Hubbard model on a two-dimensional honeycomb lattice 
by using the quantum Monte Carlo simulation with the 
stochastic series expansion (SSE). We present the results 
both in the hard-core and soft-core case. We find that 
for the hard-core case, the supersolid state is unstable to- 
wards the phase separation due to the domain wall prolif- 
eration mechanism. The transition between p = 1/2 solid 
and the superfluid is the first order one. For the soft-core 
case, due to the presence of the multiple occupation, a 
stable particle induced supersolid ( SS-p ) phase emerges 
when atoms are doped into p = 1/2 solid (i.e., fillings 
p > 1/2), while for fillings p < 1/2, the possible hole in- 
duced supersolid (SS-h ) phase is unstable towards phase 
separation. The results are qualitatively similar to the 
one obtained for a square lattice. We found the CDW-I 
to the SS-p transition is second order with the super- 



fluid density inside the SS-p scaling as p s ~ p — 1/2, 
the SS-p has the same diagonal order as the CDW-I ( 
Fig. 3). However, the SS-p to the CDW-II transition is 
first order with jumps in both p s and Pa-b- All these 
results are consistent with those achieved by the DVM 
in [ijj]. We also calculate the critical exponents of the 
transition between p = 1/2 solid and the superfluid at 
Heisenberg point. The dynamical critical exponent z and 
the correlation length exponent v we find is z = 0.15 and 
v = 0.38. Howver, w expect that z is size dependent and 
approaches zero in thermodynamic limit. This fact indi- 
cates that the solid to the SF transition is still 1st order 
even at the Heisenberg point. 

Note added. Recently, we become aware of a similar 
work [23j], where similar results are obtained. 
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